library(readr)
library(ggplot2)
library(dplyr)
library(rms)
library(synthpop)
library(tidyr)
library(plotly)
library(sjPlot)
# install the 'mi' package if you are planning on redoing imputation on the original dataset. For now the imputed dataset has already been saved
# library(mi)
source('./helper functions.R')
Reading in original dataset, performing multiple imputation and writing out imputed dataset. This chunk commented out because we’ve already done this and saved the resulting file as our starting point.
Read in dataset from Brinati paper. Dataset already been imputed based on code above.
require(readr)
brinati_covid_study_v2.imputed <- read_csv("data/brinati-covid_study_v2_imputed.csv",
col_types = cols(GENDER = col_factor(levels = c("M",
"F")), SWAB = col_factor(levels = c("0",
"1"))))
glm.orig.fit<-glm(SWAB ~ ., brinati_covid_study_v2.imputed, family='binomial')
glm.fit: fitted probabilities numerically 0 or 1 occurred
glm.orig.fit
Call: glm(formula = SWAB ~ ., family = "binomial", data = brinati_covid_study_v2.imputed)
Coefficients:
(Intercept) GENDERF AGE WBC Platelets Neutrophils Lymphocytes Monocytes Eosinophils Basophils CRP AST
-2.185e+00 -8.505e-01 1.746e-02 -1.146e-01 2.011e-03 -1.289e-01 9.374e-03 8.215e-01 -1.095e-08 -4.520e+00 5.427e-04 -7.082e-03
ALT ALP GGT LDH
1.026e-02 -1.195e-02 5.287e-03 1.001e-02
Degrees of Freedom: 278 Total (i.e. Null); 263 Residual
Null Deviance: 366.4
Residual Deviance: 239.4 AIC: 271.4
write.csv(summary(glm.orig.fit)$coef, file='output/OR-original-model.csv', row.names = FALSE)
predicted = plogis(predict(glm.orig.fit))
observed = brinati_covid_study_v2.imputed$SWAB
assessPerf(predicted, observed)
Dxy C (ROC) R2 D D:Chi-sq D:p U U:Chi-sq U:p Q Brier
7.322477e-01 8.661239e-01 5.000567e-01 4.514185e-01 1.269458e+02 NA -7.168459e-03 8.242296e-13 1.000000e+00 4.585870e-01 1.361516e-01
Intercept Slope Emax E90 Eavg S:z S:p
-2.780966e-08 1.000000e+00 6.627578e-02 5.261796e-02 2.628540e-02 -2.750989e-01 7.832403e-01
p<-plot_model(glm.orig.fit)
p<-ApplyFigureThemeLargeFontOnly(p + ylim(.1, 2.5) )
Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.
p
ggsave(filename = 'figs/OR-original-model.png', plot=p, width=6, height=6, units='in', dpi=300)

Draw ROC curve
cms.cutoffs <- lapply(seq(0.01,0.99,0.01), function(cutoff) {
ret=confusionMatrix(predicted, observed, cutoff=cutoff)
ret$cutoff = cutoff
ret
})
ff<-data.frame(t(sapply(cms.cutoffs, function(cf) c("Cutoff"=cf$cutoff, "Sensitivty"=cf$sens, "Specificity"=cf$spec))))
colnames(ff) <- c('Cutoff', 'Sensitivity','Specificity')
cutoffs_to_plot <- c(0.09, 0.3, 0.6, 0.9)
p <- ApplyFigureTheme(ff %>% mutate(fpr=1-Specificity) %>% filter(Cutoff %in% cutoffs_to_plot) %>%
ggplot(., aes(x=fpr, y=Sensitivity)) +
geom_point(size=3, color='darkblue') + xlim(0,1) + ylim(0,1)+
xlab('False Positive Rate\n(1-Specificity)') + ylab('True Positive Rate\n(Sensitivity)')
)
ggsave(filename = 'figs/example-sens-spec-on-roc.png', plot=p, width=6, height=6, units='in', dpi=300)
for (cutoff in lapply(1:length(cutoffs_to_plot), function(i) cutoffs_to_plot[1:i])) {
p <- ApplyFigureTheme(ff %>% mutate(fpr=1-Specificity) %>% filter(Cutoff %in% cutoff) %>%
ggplot(., aes(x=fpr, y=Sensitivity)) +
geom_point(size=3, color='darkblue') + xlim(0,1) + ylim(0,1)+
xlab('False Positive Rate\n(1-Specificity)') + ylab('True Positive Rate\n(Sensitivity)')
)
ggsave(filename = paste0('figs/example-sens-spec-on-roc-', cutoff[length(cutoff)], '.png'),
plot=p, width=6, height=6, units='in', dpi=300)
}
p <- ApplyFigureTheme(
ggplot(ff, aes(x=1-Specificity, y=Sensitivity)) + geom_line()+
geom_ribbon(aes(ymin = 0, ymax=Sensitivity), color=NA, fill='blue',alpha=0.3) +
xlim(0,1) + ylim(0,1)+
xlab('False Positive Rate\n(1-Specificity)') + ylab('True Positive Rate\n(Sensitivity)')
)
ggsave(filename = 'figs/example-roc.png',
plot=p, width=6, height=6, units='in', dpi=300)
p <- ff %>% mutate(fpr=1-Specificity) %>%
ggplot(., aes(x=fpr, y=Sensitivity)) + geom_point(aes(frame=Cutoff)) +
xlab('False Positive Rate (1-Specificity)') + ylab('True Positive Rate\n(Sensitivity)')
ggplotly(p)

Draw calibration
cal_breaks = seq(0, 1, .1)
cal_deciles <- data.frame(Predicted=predicted, Observed = as.integer(as.character(observed))) %>%
mutate(Bins = cut(Predicted, breaks=cal_breaks, include.lowest = TRUE))
cal_deciles_summary <- cal_deciles %>% group_by(Bins) %>% summarise(n=n(), Average_Observed = mean(Observed),
Average_Predicted=mean(Predicted))
for (cutoff in lapply(1:nrow(cal_deciles_summary), function(i) cal_deciles_summary$Bins[1:i])) {
p<-ApplyFigureThemeCalCurvePoints(
ggplot(cal_deciles_summary %>% filter(Bins %in% cutoff), aes(x=Average_Predicted, y=Average_Observed)) )
p<-ApplyFigureThemeLargeFontOnly(p)
SaveStdSquareFigure(p, filename = paste0('figs/example-cal-curve-points-', cutoff[length(cutoff)], '.png'))
}
p<-ApplyFigureThemeCalCurvePoints(ggplot(cal_deciles_summary , aes(x=Average_Predicted, y=Average_Observed)) )
p<-ApplyFigureThemeLargeFontOnly(p)
SaveStdSquareFigure(p, 'figs/cal-curve-deciles-all.png')
p<-p + geom_smooth(method='lm', se=FALSE)
SaveStdSquareFigure(p, 'figs/cal-curve-deciles-all-with-bestfit.png')
p

NA
Calculate hosmer lemeshow statistics
hl.org <- performance_hosmer(glm.orig.fit)
hl.org
# Hosmer-Lemeshow Goodness-of-Fit Test
Chi-squared: 5.333
df: 8
p-value: 0.722
Summary: model seems to fit well.
Create smoothed cal curve from Van hoorde et al
p<-CreateSmoothedCalCurvePlot(predicted, observed)
p<-ApplyFigureThemeLargeFontOnly(p)
SaveStdSquareFigure(p, 'figs/smoothed-cal-curve-orig-data.png')
p

Boot strap performance in the original dataset

Assess performance in external validation
- Look at performance as a function of prevalence
prevalences = seq(0.01, 0.99, 0.01)
brinati.syn.factory = generateSyntheticDataFactory(baseDF = brinati_covid_study_v2.imputed, method = 'boot')
dfs_prevs <- lapply(prevalences, function (prev) brinati.syn.factory(prev = prev))
dfs.prevs.perf <- lapply(dfs_prevs, function(df) {
assessPerf(predicted = plogis(predict(glm.orig.fit, newdata=as.data.frame(df))), observed = df$SWAB)
})
Calculate the different validation measures
dfs.prevs.perf.df <- data.frame(matrix(unlist(dfs.prevs.perf), nrow = length(prevalences), byrow = TRUE))
colnames(dfs.prevs.perf.df) <- names(dfs.prevs.perf[[1]])
dfs.prevs.perf.df$Prevalence = prevalences
dfs.prevs.perf.df <- gather(dfs.prevs.perf.df, Measure, Value, Dxy:`S:p`, factor_key = TRUE)
Plot the distribution of validation measures
measures = c('C (ROC)', 'Brier', 'Intercept', 'Slope')
p<-ggplot(dfs.prevs.perf.df %>% filter(Measure %in% measures), aes(x=Measure, y=Value, color=Measure)) + geom_boxplot() +
geom_point(data=data.frame(Measure = measures, Value = c(1, 0, 0, 1)), size=3, color='blue') + theme_bw()
p<-ApplyFigureThemeLargeFontOnly(p)
SaveMedRectFigure(p, 'figs/change-in-metrics-due-to-prevalence-boxplot.png')
p
p<-ggplot(dfs.prevs.perf.df %>% filter(Measure %in% measures), aes(x=Prevalence, y=Value, color=Measure, group=Measure)) + geom_line(size=1.3) + theme_bw()
p<-ApplyFigureThemeLargeFontOnly(p)
SaveMedRectFigure(p, 'figs/change-in-metrics-due-to-prevalence-linearplot.png')

p

LS0tDQp0aXRsZTogIlIgTm90ZWJvb2siDQpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sNCi0tLQ0KYGBge3IgbWVzc2FnZT1GQUxTRX0NCmxpYnJhcnkocmVhZHIpDQpsaWJyYXJ5KGdncGxvdDIpDQpsaWJyYXJ5KGRwbHlyKQ0KbGlicmFyeShybXMpDQpsaWJyYXJ5KHN5bnRocG9wKQ0KbGlicmFyeSh0aWR5cikNCmxpYnJhcnkocGxvdGx5KQ0KbGlicmFyeShzalBsb3QpDQpsaWJyYXJ5KHBlcmZvcm1hbmNlKQ0KIyBpbnN0YWxsIHRoZSAnbWknIHBhY2thZ2UgaWYgeW91IGFyZSBwbGFubmluZyBvbiByZWRvaW5nIGltcHV0YXRpb24gb24gdGhlIG9yaWdpbmFsIGRhdGFzZXQuIEZvciBub3cgdGhlIGltcHV0ZWQgZGF0YXNldCBoYXMgYWxyZWFkeSBiZWVuIHNhdmVkDQojIGxpYnJhcnkobWkpDQpzb3VyY2UoJy4vaGVscGVyIGZ1bmN0aW9ucy5SJykNCmBgYA0KDQpSZWFkaW5nIGluIG9yaWdpbmFsIGRhdGFzZXQsIHBlcmZvcm1pbmcgbXVsdGlwbGUgaW1wdXRhdGlvbiBhbmQgd3JpdGluZyBvdXQgaW1wdXRlZCBkYXRhc2V0LiBUaGlzIGNodW5rIGNvbW1lbnRlZCBvdXQgYmVjYXVzZSB3ZSd2ZSBhbHJlYWR5IGRvbmUgdGhpcyBhbmQgc2F2ZWQgdGhlIHJlc3VsdGluZyBmaWxlIGFzIG91ciBzdGFydGluZyBwb2ludC4NCmBgYHtyIGVjaG89RkFMU0UsIHdhcm5pbmc9RkFMU0V9DQojIA0KIyBicmluYXRpX2NvdmlkX3N0dWR5X3YyIDwtIHJlYWRfY3N2KCJkYXRhL2JyaW5hdGktY292aWRfc3R1ZHlfdjIuY3N2IiwNCiMgY29sX3R5cGVzID0gY29scyhHRU5ERVIgPSBjb2xfZmFjdG9yKGxldmVscyA9IGMoIk0iLA0KIyAiRiIpKSwgU1dBQiA9IGNvbF9mYWN0b3IobGV2ZWxzID0gYygiMCIsDQojICIxIikpKSkNCiMgYnJpbmF0aV9jb3ZpZF9zdHVkeV92MiA8LSBicmluYXRpX2NvdmlkX3N0dWR5X3YyICU+JSBtdXRhdGVfaWYoaXMubnVtZXJpYywgLmZ1bnMgPSBmdW5jdGlvbih4KSB7eCswLjAwMX0pDQojIGJyYWluYXRpLm1pIDwtIG1pc3NpbmdfZGF0YS5mcmFtZShhcy5kYXRhLmZyYW1lKGJyaW5hdGlfY292aWRfc3R1ZHlfdjIpLCBmYXZvcl9wb3NpdGl2ZSA9IFRSVUUpDQojIA0KIyANCiMgYnJhaW5hdGkubWk8LWNoYW5nZShicmFpbmF0aS5taSwgeT1jKCdMeW1waG9jeXRlcycsICdCYXNvcGhpbHMnLCAnTW9ub2N5dGVzJywnRW9zaW5vcGhpbHMnLCAnQmFzb3BoaWxzJyksIHdoYXQ9J3R5cGUnLCB0bz0ncG9zaXRpdmUtY29udGludW91cycpDQojIHNob3coYnJhaW5hdGkubWkpDQojIGltYWdlKGJyYWluYXRpLm1pKQ0KIyANCiMgb3B0aW9ucyhtYy5jb3JlcyA9IDQpDQojIGltcHV0YXRpb25zIDwtIG1pKGJyYWluYXRpLm1pLCBuLml0ZXIgPSA5MCwgbi5jaGFpbnMgPSA0LCBtYXgubWludXRlcyA9IDIwKSANCiMgc2hvdyhpbXB1dGF0aW9ucykNCiMgcm91bmQobWlwcGx5KGltcHV0YXRpb25zLCBtZWFuLCB0by5tYXRyaXggPSBUUlVFKSwgMykNCiMgUmhhdHMoaW1wdXRhdGlvbnMpDQojIGRmcyA8LSBjb21wbGV0ZShpbXB1dGF0aW9ucywgbT0yKQ0KIyBicmluYXRpX2NvdmlkX3YyLmltcHV0ZWQgPC0gZGZzW1sxXV1bLDE6MTZdDQojIHdyaXRlLmNzdihicmluYXRpX2NvdmlkX3YyLmltcHV0ZWQsIGZpbGU9Jy4vZGF0YS9icmluYXRpLWNvdmlkX3N0dWR5X3YyX2ltcHV0ZWQuY3N2Jywgcm93Lm5hbWVzID0gRkFMU0UpDQpgYGANCg0KUmVhZCBpbiBkYXRhc2V0IGZyb20gQnJpbmF0aSBwYXBlci4gRGF0YXNldCBhbHJlYWR5IGJlZW4gaW1wdXRlZCBiYXNlZCBvbiBjb2RlIGFib3ZlLiANCmBgYHtyfQ0KcmVxdWlyZShyZWFkcikNCmJyaW5hdGlfY292aWRfc3R1ZHlfdjIuaW1wdXRlZCA8LSByZWFkX2NzdigiZGF0YS9icmluYXRpLWNvdmlkX3N0dWR5X3YyX2ltcHV0ZWQuY3N2IiwNCiAgY29sX3R5cGVzID0gY29scyhHRU5ERVIgPSBjb2xfZmFjdG9yKGxldmVscyA9IGMoIk0iLA0KICAiRiIpKSwgU1dBQiA9IGNvbF9mYWN0b3IobGV2ZWxzID0gYygiMCIsDQogICIxIikpKSkNCmBgYA0KDQoNCmBgYHtyfQ0KDQpnbG0ub3JpZy5maXQ8LWdsbShTV0FCIH4gLiwgYnJpbmF0aV9jb3ZpZF9zdHVkeV92Mi5pbXB1dGVkLCBmYW1pbHk9J2Jpbm9taWFsJykNCmdsbS5vcmlnLmZpdA0Kd3JpdGUuY3N2KHN1bW1hcnkoZ2xtLm9yaWcuZml0KSRjb2VmLCBmaWxlPSdvdXRwdXQvT1Itb3JpZ2luYWwtbW9kZWwuY3N2Jywgcm93Lm5hbWVzID0gRkFMU0UpDQpwcmVkaWN0ZWQgPSBwbG9naXMocHJlZGljdChnbG0ub3JpZy5maXQpKQ0Kb2JzZXJ2ZWQgPSBicmluYXRpX2NvdmlkX3N0dWR5X3YyLmltcHV0ZWQkU1dBQg0KYXNzZXNzUGVyZihwcmVkaWN0ZWQsIG9ic2VydmVkKQ0KcDwtcGxvdF9tb2RlbChnbG0ub3JpZy5maXQpDQpwPC1BcHBseUZpZ3VyZVRoZW1lTGFyZ2VGb250T25seShwICsgeWxpbSguMSwgMi41KSApIA0KICANCnANCmdnc2F2ZShmaWxlbmFtZSA9ICdmaWdzL09SLW9yaWdpbmFsLW1vZGVsLnBuZycsIHBsb3Q9cCwgd2lkdGg9NiwgaGVpZ2h0PTYsIHVuaXRzPSdpbicsIGRwaT0zMDApDQpgYGANCg0KIyBEcmF3IFJPQyBjdXJ2ZQ0KYGBge3J9DQpjbXMuY3V0b2ZmcyA8LSBsYXBwbHkoc2VxKDAuMDEsMC45OSwwLjAxKSwgZnVuY3Rpb24oY3V0b2ZmKSB7DQogIHJldD1jb25mdXNpb25NYXRyaXgocHJlZGljdGVkLCBvYnNlcnZlZCwgY3V0b2ZmPWN1dG9mZikNCiAgcmV0JGN1dG9mZiA9IGN1dG9mZg0KICByZXQNCiAgfSkNCmZmPC1kYXRhLmZyYW1lKHQoc2FwcGx5KGNtcy5jdXRvZmZzLCBmdW5jdGlvbihjZikgYygiQ3V0b2ZmIj1jZiRjdXRvZmYsICJTZW5zaXRpdnR5Ij1jZiRzZW5zLCAiU3BlY2lmaWNpdHkiPWNmJHNwZWMpKSkpDQpjb2xuYW1lcyhmZikgPC0gYygnQ3V0b2ZmJywgJ1NlbnNpdGl2aXR5JywnU3BlY2lmaWNpdHknKQ0KDQoNCmN1dG9mZnNfdG9fcGxvdCA8LSBjKDAuMDksIDAuMywgMC42LCAwLjkpDQpwIDwtIEFwcGx5RmlndXJlVGhlbWUoZmYgJT4lIG11dGF0ZShmcHI9MS1TcGVjaWZpY2l0eSkgJT4lIGZpbHRlcihDdXRvZmYgJWluJSBjdXRvZmZzX3RvX3Bsb3QpICU+JQ0KICBnZ3Bsb3QoLiwgYWVzKHg9ZnByLCB5PVNlbnNpdGl2aXR5KSkgKyANCiAgICBnZW9tX3BvaW50KHNpemU9MywgY29sb3I9J2RhcmtibHVlJykgKyB4bGltKDAsMSkgKyB5bGltKDAsMSkrDQogICAgeGxhYignRmFsc2UgUG9zaXRpdmUgUmF0ZVxuKDEtU3BlY2lmaWNpdHkpJykgKyB5bGFiKCdUcnVlIFBvc2l0aXZlIFJhdGVcbihTZW5zaXRpdml0eSknKQ0KICApDQpnZ3NhdmUoZmlsZW5hbWUgPSAnZmlncy9leGFtcGxlLXNlbnMtc3BlYy1vbi1yb2MucG5nJywgcGxvdD1wLCB3aWR0aD02LCBoZWlnaHQ9NiwgdW5pdHM9J2luJywgZHBpPTMwMCkNCg0KZm9yIChjdXRvZmYgaW4gbGFwcGx5KDE6bGVuZ3RoKGN1dG9mZnNfdG9fcGxvdCksIGZ1bmN0aW9uKGkpIGN1dG9mZnNfdG9fcGxvdFsxOmldKSkgew0KICBwIDwtIEFwcGx5RmlndXJlVGhlbWUoZmYgJT4lIG11dGF0ZShmcHI9MS1TcGVjaWZpY2l0eSkgJT4lIGZpbHRlcihDdXRvZmYgJWluJSBjdXRvZmYpICU+JQ0KICAgIGdncGxvdCguLCBhZXMoeD1mcHIsIHk9U2Vuc2l0aXZpdHkpKSArIA0KICAgICAgZ2VvbV9wb2ludChzaXplPTMsIGNvbG9yPSdkYXJrYmx1ZScpICsgeGxpbSgwLDEpICsgeWxpbSgwLDEpKw0KICAgICAgeGxhYignRmFsc2UgUG9zaXRpdmUgUmF0ZVxuKDEtU3BlY2lmaWNpdHkpJykgKyB5bGFiKCdUcnVlIFBvc2l0aXZlIFJhdGVcbihTZW5zaXRpdml0eSknKQ0KICApDQogIGdnc2F2ZShmaWxlbmFtZSA9IHBhc3RlMCgnZmlncy9leGFtcGxlLXNlbnMtc3BlYy1vbi1yb2MtJywgY3V0b2ZmW2xlbmd0aChjdXRvZmYpXSwgICcucG5nJyksIA0KICAgICAgICAgcGxvdD1wLCB3aWR0aD02LCBoZWlnaHQ9NiwgdW5pdHM9J2luJywgZHBpPTMwMCkNCn0NCnAgPC0gQXBwbHlGaWd1cmVUaGVtZSgNCiAgZ2dwbG90KGZmLCBhZXMoeD0xLVNwZWNpZmljaXR5LCB5PVNlbnNpdGl2aXR5KSkgKyBnZW9tX2xpbmUoKSsgDQogICAgZ2VvbV9yaWJib24oYWVzKHltaW4gPSAwLCB5bWF4PVNlbnNpdGl2aXR5KSwgY29sb3I9TkEsIGZpbGw9J2JsdWUnLGFscGhhPTAuMykgKyANCiAgICB4bGltKDAsMSkgKyB5bGltKDAsMSkrDQogICAgICB4bGFiKCdGYWxzZSBQb3NpdGl2ZSBSYXRlXG4oMS1TcGVjaWZpY2l0eSknKSArIHlsYWIoJ1RydWUgUG9zaXRpdmUgUmF0ZVxuKFNlbnNpdGl2aXR5KScpIA0KKQ0KZ2dzYXZlKGZpbGVuYW1lID0gJ2ZpZ3MvZXhhbXBsZS1yb2MucG5nJywgDQogICAgICAgICBwbG90PXAsIHdpZHRoPTYsIGhlaWdodD02LCB1bml0cz0naW4nLCBkcGk9MzAwKQ0KcCA8LSBmZiAlPiUgbXV0YXRlKGZwcj0xLVNwZWNpZmljaXR5KSAlPiUNCiAgZ2dwbG90KC4sIGFlcyh4PWZwciwgeT1TZW5zaXRpdml0eSkpICsgZ2VvbV9wb2ludChhZXMoZnJhbWU9Q3V0b2ZmKSkgKw0KICAgIHhsYWIoJ0ZhbHNlIFBvc2l0aXZlIFJhdGUgKDEtU3BlY2lmaWNpdHkpJykgKyB5bGFiKCdUcnVlIFBvc2l0aXZlIFJhdGVcbihTZW5zaXRpdml0eSknKQ0KICANCiAgDQpnZ3Bsb3RseShwKQ0KYGBgDQoNCiMgRHJhdyBjYWxpYnJhdGlvbg0KYGBge3J9DQpjYWxfYnJlYWtzID0gc2VxKDAsIDEsIC4xKQ0KY2FsX2RlY2lsZXMgPC0gZGF0YS5mcmFtZShQcmVkaWN0ZWQ9cHJlZGljdGVkLCBPYnNlcnZlZCA9IGFzLmludGVnZXIoYXMuY2hhcmFjdGVyKG9ic2VydmVkKSkpICU+JQ0KICBtdXRhdGUoQmlucyA9IGN1dChQcmVkaWN0ZWQsIGJyZWFrcz1jYWxfYnJlYWtzLCBpbmNsdWRlLmxvd2VzdCA9IFRSVUUpKQ0KY2FsX2RlY2lsZXNfc3VtbWFyeSA8LSBjYWxfZGVjaWxlcyAlPiUgZ3JvdXBfYnkoQmlucykgJT4lIHN1bW1hcmlzZShuPW4oKSwgQXZlcmFnZV9PYnNlcnZlZCA9IG1lYW4oT2JzZXJ2ZWQpLA0KICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBBdmVyYWdlX1ByZWRpY3RlZD1tZWFuKFByZWRpY3RlZCkpIA0KDQoNCmZvciAoY3V0b2ZmIGluIGxhcHBseSgxOm5yb3coY2FsX2RlY2lsZXNfc3VtbWFyeSksIGZ1bmN0aW9uKGkpIGNhbF9kZWNpbGVzX3N1bW1hcnkkQmluc1sxOmldKSkgew0KICBwPC1BcHBseUZpZ3VyZVRoZW1lQ2FsQ3VydmVQb2ludHMoDQogICAgZ2dwbG90KGNhbF9kZWNpbGVzX3N1bW1hcnkgJT4lIGZpbHRlcihCaW5zICVpbiUgY3V0b2ZmKSwgYWVzKHg9QXZlcmFnZV9QcmVkaWN0ZWQsIHk9QXZlcmFnZV9PYnNlcnZlZCkpICkNCiAgcDwtQXBwbHlGaWd1cmVUaGVtZUxhcmdlRm9udE9ubHkocCkNCiAgU2F2ZVN0ZFNxdWFyZUZpZ3VyZShwLCBmaWxlbmFtZSA9IHBhc3RlMCgnZmlncy9leGFtcGxlLWNhbC1jdXJ2ZS1wb2ludHMtJywgY3V0b2ZmW2xlbmd0aChjdXRvZmYpXSwgICcucG5nJykpDQp9DQoNCnA8LUFwcGx5RmlndXJlVGhlbWVDYWxDdXJ2ZVBvaW50cyhnZ3Bsb3QoY2FsX2RlY2lsZXNfc3VtbWFyeSAsIGFlcyh4PUF2ZXJhZ2VfUHJlZGljdGVkLCB5PUF2ZXJhZ2VfT2JzZXJ2ZWQpKSApDQpwPC1BcHBseUZpZ3VyZVRoZW1lTGFyZ2VGb250T25seShwKQ0KDQpTYXZlU3RkU3F1YXJlRmlndXJlKHAsICdmaWdzL2NhbC1jdXJ2ZS1kZWNpbGVzLWFsbC5wbmcnKQ0KDQpwPC1wICsgZ2VvbV9zbW9vdGgobWV0aG9kPSdsbScsIHNlPUZBTFNFKQ0KU2F2ZVN0ZFNxdWFyZUZpZ3VyZShwLCAnZmlncy9jYWwtY3VydmUtZGVjaWxlcy1hbGwtd2l0aC1iZXN0Zml0LnBuZycpDQpwDQogICAgICAgICAgICANCmBgYA0KDQpDYWxjdWxhdGUgaG9zbWVyIGxlbWVzaG93IHN0YXRpc3RpY3MNCmBgYHtyfQ0KDQpobC5vcmcgPC0gcGVyZm9ybWFuY2VfaG9zbWVyKGdsbS5vcmlnLmZpdCkNCmhsLm9yZw0KYGBgDQoNCkNyZWF0ZSBzbW9vdGhlZCBjYWwgY3VydmUgZnJvbSBWYW4gaG9vcmRlIGV0IGFsDQpgYGB7cn0NCnA8LUNyZWF0ZVNtb290aGVkQ2FsQ3VydmVQbG90KHByZWRpY3RlZCwgb2JzZXJ2ZWQpDQpwPC1BcHBseUZpZ3VyZVRoZW1lTGFyZ2VGb250T25seShwKQ0KU2F2ZVN0ZFNxdWFyZUZpZ3VyZShwLCAnZmlncy9zbW9vdGhlZC1jYWwtY3VydmUtb3JpZy1kYXRhLnBuZycpDQpwDQpgYGANCg0KQm9vdCBzdHJhcCBwZXJmb3JtYW5jZSBpbiB0aGUgb3JpZ2luYWwgZGF0YXNldA0KYGBge3Igd2FybmluZz1GQUxTRSwgZWNobz1GQUxTRX0NCmRmIDwtIGFzLmRhdGEuZnJhbWUoYnJpbmF0aV9jb3ZpZF9zdHVkeV92Mi5pbXB1dGVkKQ0KbkJvb3QgPSAxMDANCnNldC5zZWVkKDEzNDIzNDkpDQp0cmFpbklkeEJvb3QgPC0gbGFwcGx5KDE6bkJvb3QsIGZ1bmN0aW9uKGkpIHNhbXBsZSgxOm5yb3coZGYpLCBzaXplPW5yb3coZGYpLCByZXBsYWNlPVRSVUUpKQ0KDQpib290UGVyZi5vcmlnZGF0YSA8LSBsYXBwbHkodHJhaW5JZHhCb290LCANCiAgICAgICAgICAgICAgICAgICBmdW5jdGlvbih0cmFpbklkeCkgDQogICAgICAgICAgICAgICAgICAgICBhc3Nlc3NTaW5nbGVUcmFpblRlc3QodHJhaW5JZHgsICgxOm5yb3coZGYpKVstdHJhaW5JZHhdLCBkZiwgZ2xtLm1vZGVsID0gU1dBQn4gLiwgb3V0Y29tZS52YXIubmFtZSA9ICdTV0FCJykNCiAgICAgICAgICAgICAgICAgICApIA0KDQptZWFzdXJlcyA9IGMoJ0MgKFJPQyknLCAnQnJpZXInLCAnSW50ZXJjZXB0JywgJ1Nsb3BlJykNCmJyaW5hdGkub3JpZy5ib290LnBlcmYgPC0gZGF0YS5mcmFtZShtYXRyaXgoDQogIHVubGlzdChsYXBwbHkoYm9vdFBlcmYub3JpZ2RhdGEsIGZ1bmN0aW9uKGJvb3QpIGJvb3QkdGVzdC5wZXJmKSksIA0KICBucm93ID0gbkJvb3QsIGJ5cm93ID0gVFJVRSkpDQpjb2xuYW1lcyhicmluYXRpLm9yaWcuYm9vdC5wZXJmKSA8LSBuYW1lcyhib290UGVyZi5vcmlnZGF0YVtbMV1dJHRlc3QucGVyZikNCg0KYnJpbmF0aS5vcmlnLmJvb3QucGVyZiA8LSBnYXRoZXIoYnJpbmF0aS5vcmlnLmJvb3QucGVyZiwgTWVhc3VyZSwgVmFsdWUsIER4eTpgUzpwYCwgZmFjdG9yX2tleSA9IFRSVUUpDQpwPC1nZ3Bsb3QoYnJpbmF0aS5vcmlnLmJvb3QucGVyZiAlPiUgZmlsdGVyKE1lYXN1cmUgJWluJSBtZWFzdXJlcyksIGFlcyh4PU1lYXN1cmUsIHk9VmFsdWUsICBjb2xvcj1NZWFzdXJlKSkgKyBnZW9tX2JveHBsb3QoKSArDQogIGdlb21fcG9pbnQoZGF0YT1kYXRhLmZyYW1lKE1lYXN1cmUgPSBtZWFzdXJlcywgVmFsdWUgPSBjKDEsIDAsIDAsIDEpKSwgc2l6ZT0zLCBjb2xvcj0nYmx1ZScpICsgdGhlbWVfYncoKQ0KcDwtQXBwbHlGaWd1cmVUaGVtZShwKQ0KU2F2ZVN0ZFNxdWFyZUZpZ3VyZShwLCAnZmlncy9vcmlnLW1vZGVsLWJvb3RzdHJhcC1ldmFsdWF0aW9uLnBuZycpDQpwDQpgYGANCg0KIyBBc3Nlc3MgcGVyZm9ybWFuY2UgaW4gZXh0ZXJuYWwgdmFsaWRhdGlvbg0KDQpJSS4gTG9vayBhdCBwZXJmb3JtYW5jZSBhcyBhIGZ1bmN0aW9uIG9mIHByZXZhbGVuY2UNCmBgYHtyfQ0KcHJldmFsZW5jZXMgPSBzZXEoMC4wMSwgMC45OSwgMC4wMSkNCmJyaW5hdGkuc3luLmZhY3RvcnkgPSBnZW5lcmF0ZVN5bnRoZXRpY0RhdGFGYWN0b3J5KGJhc2VERiA9IGJyaW5hdGlfY292aWRfc3R1ZHlfdjIuaW1wdXRlZCwgbWV0aG9kID0gJ2Jvb3QnKQ0KZGZzX3ByZXZzIDwtIGxhcHBseShwcmV2YWxlbmNlcywgZnVuY3Rpb24gKHByZXYpIGJyaW5hdGkuc3luLmZhY3RvcnkocHJldiA9IHByZXYpKQ0KZGZzLnByZXZzLnBlcmYgPC0gbGFwcGx5KGRmc19wcmV2cywgZnVuY3Rpb24oZGYpIHsNCiAgYXNzZXNzUGVyZihwcmVkaWN0ZWQgPSBwbG9naXMocHJlZGljdChnbG0ub3JpZy5maXQsIG5ld2RhdGE9YXMuZGF0YS5mcmFtZShkZikpKSwgb2JzZXJ2ZWQgPSBkZiRTV0FCKQ0KfSkNCmBgYA0KDQpDYWxjdWxhdGUgdGhlIGRpZmZlcmVudCB2YWxpZGF0aW9uIG1lYXN1cmVzDQpgYGB7cn0NCg0KZGZzLnByZXZzLnBlcmYuZGYgPC0gZGF0YS5mcmFtZShtYXRyaXgodW5saXN0KGRmcy5wcmV2cy5wZXJmKSwgbnJvdyA9IGxlbmd0aChwcmV2YWxlbmNlcyksIGJ5cm93ID0gVFJVRSkpDQpjb2xuYW1lcyhkZnMucHJldnMucGVyZi5kZikgPC0gbmFtZXMoZGZzLnByZXZzLnBlcmZbWzFdXSkNCmRmcy5wcmV2cy5wZXJmLmRmJFByZXZhbGVuY2UgPSBwcmV2YWxlbmNlcw0KZGZzLnByZXZzLnBlcmYuZGYgPC0gZ2F0aGVyKGRmcy5wcmV2cy5wZXJmLmRmLCBNZWFzdXJlLCBWYWx1ZSwgRHh5OmBTOnBgLCBmYWN0b3Jfa2V5ID0gVFJVRSkNCmBgYA0KDQpQbG90IHRoZSBkaXN0cmlidXRpb24gb2YgdmFsaWRhdGlvbiBtZWFzdXJlcw0KYGBge3J9DQptZWFzdXJlcyA9IGMoJ0MgKFJPQyknLCAnQnJpZXInLCAnSW50ZXJjZXB0JywgJ1Nsb3BlJykNCnA8LWdncGxvdChkZnMucHJldnMucGVyZi5kZiAlPiUgZmlsdGVyKE1lYXN1cmUgJWluJSBtZWFzdXJlcyksIGFlcyh4PU1lYXN1cmUsIHk9VmFsdWUsICBjb2xvcj1NZWFzdXJlKSkgKyBnZW9tX2JveHBsb3QoKSArDQogIGdlb21fcG9pbnQoZGF0YT1kYXRhLmZyYW1lKE1lYXN1cmUgPSBtZWFzdXJlcywgVmFsdWUgPSBjKDEsIDAsIDAsIDEpKSwgc2l6ZT0zLCBjb2xvcj0nYmx1ZScpICsgdGhlbWVfYncoKQ0KcDwtQXBwbHlGaWd1cmVUaGVtZUxhcmdlRm9udE9ubHkocCkNClNhdmVNZWRSZWN0RmlndXJlKHAsICdmaWdzL2NoYW5nZS1pbi1tZXRyaWNzLWR1ZS10by1wcmV2YWxlbmNlLWJveHBsb3QucG5nJykNCnANCnA8LWdncGxvdChkZnMucHJldnMucGVyZi5kZiAlPiUgZmlsdGVyKE1lYXN1cmUgJWluJSBtZWFzdXJlcyksIGFlcyh4PVByZXZhbGVuY2UsIHk9VmFsdWUsICBjb2xvcj1NZWFzdXJlLCBncm91cD1NZWFzdXJlKSkgKyBnZW9tX2xpbmUoc2l6ZT0xLjMpICsgdGhlbWVfYncoKQ0KcDwtQXBwbHlGaWd1cmVUaGVtZUxhcmdlRm9udE9ubHkocCkNClNhdmVNZWRSZWN0RmlndXJlKHAsICdmaWdzL2NoYW5nZS1pbi1tZXRyaWNzLWR1ZS10by1wcmV2YWxlbmNlLWxpbmVhcnBsb3QucG5nJykNCnANCmBgYA0KDQo=